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Abstract. The essential features of a full potential electronic structure method 
using Linear Muffin- Tin Orbitals (LMTOs) are presented. The electron density 
and potential in the this method are represented with no inherent geometrical 
approximation. This method allows the calculation of total energies and forces with 
arbitrary accuracy while sacrificing much of the efficiency and physical content of 
approximate methods such as the LMTO-ASA method. 



1 Introduction 

This paper describes a particular implementation of a full-potential electronic 
structure method using Linear Muffin-Tin Orbitals (LMTO's) [[TP] as basis 
functions. There have been several "FP-LMTO" implementations||,U|,|,0| . 
The one described here has not been published in detail, although calculations 
performed with this method have been reported for quite some time. There 
are many aspects to an electronic structure method. This paper is focussed 
on those aspects which enable a full potential treatment. Relatively small 
details pertaining to full-potential methods will be discussed while larger 
details having to do with, for example, relativity will not be. 

The emphasis of a variational full-potential method is somewhat different 
from that of a method such as the LMTO-ASA method. The emphasis of 
the former is on the completeness of the basis while in the latter it is in the 
physical content (and interpretability) of the basis. These concepts are, of 
course, intimately related, but the emphasis is different. 

The exposition here is for an infinite system periodic in three dimensions. 
This method has been implemented for two-dimensional systems, Q but that 
will not be discussed here. 



Notation 

Papers on electronic structure methods unavoidably carry a high overhead 
in functional symbols and indices. It is simplest to define here, without moti- 
vation, the special symbols and functions that will be used in this paper, for 
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future reference. These special functions (although not necessarily the sym- 
bols used here) have been used extensively in LMTO documentation and are 
largely due to Andersen. Qj 



Spherical harmonics: 

y tm (r) = \ l Y lm {r) (1) 

Qm(f) = {mi y ' m(f) (2) 

C lm (r) = i e C lm (f) (3) 
where Y is a spherical harmonic. |j| 

Bessel functions: 

K t {K,r) = -K^{ ni{ « r) -' lk{Rr) i <0 (4) 
I ni(Hir) k 2 > 

K L (K,r)=Ki(K,r)y L (t) (5) 
Jt(K,r) = j e (nr)/n e (6) 
J L (n,r) =Jt[K,r)y L {r) (7) 



where L denotes Im and and jg are spherical Neumann and Bessel func- 
tions, respectively. 



Geometry: For computational purposes, the crystal is divided into non- 
overlapping spheres surrounding atomic sites (muffin-tin spheres) where the 
charge density and potential vary rapidly and the interstitial region between 
the spheres, where the charge density and potential vary slowly. This is the 
muffin-tin geometry used as an idealized potential and charge density in 
early electronic structure methods (KKR and APW). Here, the division is a 
computational one, and does not restrict the final shape of the charge density 
or potential. In the muffin-tin spheres, the basis functions, electron density, 
and potential are expanded in spherical waves; in the interstitial region, the 
basis functions, electron density, and potential are expanded in Fourier series. 

There are many relevant considerations in choosing muffin-tin radii. As- 
suming all expansions are taken to convergence, the density and potential 
depend on the muffin-tin radii only through the dependence of basis func- 
tions on the radii. As discussed below, basis functions have a different func- 
tional form inside the muffin-tin spheres, and the choice of muffin-tin radius 
affects this crossover. Hence, assuming the Hamiltonian is the same inside 
and outside the spheres (the treatment of relativity may affect this as dis- 
cussed below), the muffin-tin radii are variational parameters and the opti- 
mum choice minimizes the total energy. If the basis is large enough however 
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(suitably complete within and without the spheres), the energy is insensitive 
to the choice of radii. A reasonable choice results from choosing radii that are 
within both the minimum in charge density and the maximum in potential 
along a line between nearest neighbors. Relativistic effects are usually taken 
into account only in the muffin-tin spheres, in which case the Hamiltonian 
depends on the radii; hence when relativistic effects are important, the radii 
are not variational parameters. 

In what follows, lattice positions are vectors R — Rn, integer multiples 
of a basis R. Atomic positions in the unit cell are denoted by r. A set of 
atomic positions invariant under the point group of the lattice are said to be 
of the same symmetry type, t. Similarly, in the reciprocal lattice, vectors are 
g = Gn for the reciprocal basis G = 2irR~ T . Brillouin zone (or reciprocal 
unit cell) vectors are denoted by k. 



Symmetric functions: Within the muffin-tin region, functions invariant 
are expressed in harmonic series. If f(r) is such a function, at site r 



fir) =Y,fht{r T )D ht {V T r T ) 
Dht(r) = *y]aht(rn)Ct hm (r) 



(8a) 
(8b) 



In Equation (pa]), T> T is a transformation to a coordinate system local to site 
r; the local coordinates of sites of the same type are related by an element 
of the crystal point group that takes one site into another. Expressed in this 
way, the functional form of Dht (Equation (|Sb|)) depends only on symmetry 
type. 

In the interstitial region, symmetric functions are expressed in Fourier 
series: 



fir) 



Y,f(S)D s (r) 



D s {r) = Y, elgr 

965 



(9a) 
(9b) 



The sum in Equation (pa|) is over symmetry stars S of the reciprocal lattice. 



2 Basis Set 



2.1 Interstitial 



In the interstitial region (symbolically X) between the muffin-tin spheres, 
bases are Bloch sums of spherical Hankel or Neumann functions: 



rex 



J2e lk - R IC u ( Kl , Ir-Ti-RDy^iVrtir-Ti-R)) (10) 
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The rotation T> T in ( |l0| ) takes the argument into a coordinate system local 
to each site r. The purpose of this will be made evident later. The function 
on the right hand side of Equation ( |l0| ) is sometimes called the envelope 
function. 

Notice the parameters, specifying a basis function, inherent in this def- 
inition. They are the site r in the unit cell on which the spherical wave is 
based, the angular momentum parameters I and m of the spherical wave 
with respect to its parent cell, and the kinetic energy n 2 of the basis in the 
interstitial region. The angular momentum parameters specifying the basis 
set are chosen to represent the atomic states from which crystal eigenstates 
are derived. In the LMTO-ASA, it is usual to include £ bases one higher 
than the highest relevant band. In the method described here, this is rarely 
necessary, possibly because of the multiplicity of bases with the same angular 
momentum parameters. It is usual to use "multiple k" basis sets, having all 
parameters except the tail parameter the same. 

There appears to be no simple algorithm for choosing a good set of in- 
terstitial kinetic energy parameters. Schemes such as bracketing the relevant 
energy spectrum have been proposed. The optimum set would minimize 
the total energy. This can be done but is time consuming even for relatively 
simple systems. It seems, however that parameter sets obtained in this way for 
simple systems in representative configurations can give good results when 
used for related systems over a broad pressure range. Thus good sets are 
arrived at through some experimentation. The choice can be important as 
it's possible to pick a set of parameters that will give very bad results, and 
the parameter set used in any new calculation should be always checked for 
stability. 

2.2 Muffin Tins 

In the muffin-tin spheres, bases are linear combinations of spherical waves 
matching continuously and differentiably to the envelope function at the 
muffin-tin sphere. The envelope function IC may be expanded in a series of 
spherical Bessel functions about any site except it's center. A basis function 
on a muffin-tin sphere in the unit cell at R = is therefore 

^(fc,r)| =Y / e tkR Y / yL(^rr T )( s T )S(R, 0)S(t, n)S(L, L t ) 

+ jL(K,S T )B LtLi (Ki,T-T'-R 

= y^rr) ( Ki{Ki,a T )S(r, n)S(L, U) 

L 

+ Jh{n,s r )B LiLi {K, t ,T-T' ',k) 



(11) 
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where r T = r — r and _B is equivalent to the KKR structure constant. [[10[ 
The unitary transformation applied to B rotates components into site-local 
coordinates from the left and right. 

Equation ( pd| ) is compactly expressed by defining a two-component row 
vector K so that 

K t {K,r) = (K i (K,r),Jt{K,r)) (12) 
and a two component column vector S so that 

s «'<— • *> - U£££#,) ■ < i3 > 

Then the value of a basis function on a muffin-tin boundary is expressed 
simply as 

&(fc,r) = V3'L(Prrr)K/(K i ,s T )SL.L i (/e il T-T',fe) (14) 

r T =s T — 

The radial part a basis function inside a muffin-tin sphere is a linear 
combination of atomic like functions <fi and their energy derivatives (j) fl[H 
matching continuously and differentiably to the radial function K in Equation 
14[) . Collecting <f> and ^ in a row vector 



U(e,r) = (<£(e,r),0(e,r)) 



a simple case of this matching condition may be expressed as U(e, s)fl(e, k) — 
K(k, s) and U'(e, s)Q(e, k) = K'(k, s), where Q is a matrix of order 2. 

The use of these radial functions in the method described here is different 
than that used by most other methods, however. For the broadest utility, 
a basis set must be flexible enough to describe energy levels derived from 
atomic states having different principle quantum numbers but the same an- 
gular momentum quantum number. For example, describing the properties 
of elemental actinides at any pressure requires a basis with both 6p and 7p 
character. Similarly, an adequate calculation of the structural properties of 
transition metal oxides requires both semi-core and valence s and p states on 
the transition metal ions. The description of the evolution of core states from 
localized to itinerant under pressure also requires multiple principle quantum 
numbers per I value. It is usual in LMTO-based methods to perform calcu- 
lations for the eigenstates and eigenvalues of "semi-core" and valence states 
separately, using a different basis set, with a single set of energy parameters 
{e^}, for each "energy panel". This approach fails when energy panels over- 
lap, and has the disadvantage that the set of eigenvectors is not an orthogonal 
set. The problem of "ghost bands" also arises. || 

In the method described here, bases corresponding to multiple princi- 
ple quantum numbers are contained within a single, fully hybridizing ba- 
sis set. This is accomplished simply by using functions <$> and <f> calculated 
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with energies {e n i} corresponding to different principal quantum numbers 
n to describe the radial dependence of a basis in the muffin-tin spheres. 
The Hamiltonian matrix for an actinide, for example, will have elements 
(^6p|^|V'7p) an d the overlap matrix elements (V'6p|V'7p)j We may formally 
express the radial part of basis i in a muffin-tin sphere by the function 
/(r) = ai(n£)\J(e n i , s)Q(e n i, Ki) but in practice it is sufficient to restrict 
the coefficients by a,i{n£) = 8(n, rii) so that the basis set (although not eigen- 
vectors) will have pure principal quantum number "parentage" . This method 
of expanding the energy range of a basis set has been used (and reported) 
extensively. Representative calculations in which this method was essential 
are described in Reference Jll| . 

Thus another parameter specifying a basis function is the set of energy 
parameters {eu} that will be used to calculate the radial basis functions 
4>tz and 4>te used to express the basis function in muffin-tin spheres of each 
symmetry type. A basis function in a muffin-tin sphere is therefore 

^i{k,r) = V* UtL(ei,V T r T )f2u(ei,K,i)SL,L i (KuT-T , ,k) (15) 

r T <s t *-T ' 

where means "use the energy parameter e n i corresponding to the principal 
quantum number specified for basis i" and 

U ti (e,r)=y£(f)Utf(e,r) . (16) 

The necessary cutoff in angular momentum has now been made explicit. The 
2x2 matrix J? matches U to K continuously and differentiably at the muffin- 
tin radius. Specifically, Q is specified by 

Me, s t ) Me, s t )\ ( = / Ke(K, s t ) », s t )\ 
M^ st) Me, ^) J n,) \lC' £ (K,s t ) Jl(K,s t )J v > 

In principle, and as programmed, each (t£k) basis can use its own unique 
energy set. It is more usual to use a common energy set for a set of basis 
states giving rise to bands of similar energy within the scope of a particular 
calculation. The configuration of the basis shown in Table [l] for example uses 
a set of energies for "semi-core" 6s and 6p bases, and another set of energies to 
represent "valence" bases. The calculation of energies in an energy parameter 
set is discussed below. 

A parameter introduced in (|l^) is the angular momentum cut-off £ m . In 
most cases, a converged total energy is achieved with values £ m ~ 6 — 8. Note 
that since a basis set generally contains functions based on spherical waves 
with £ < 3, the KKR structure constant in (Oh is rectangular. 



3 Matrix Elements 
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Table 1. Parameters for typical basis set for an elemental actinide: parent angular 
momentum parameter (£), energy set for radial expansions (e-set), and the index 
of the kinetic energy in the interstitial region(K-index). A typical set of k 2 values, 
corresponding to the kinetic energy indices, is given at the bottom of the table. 
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-1.96582916 
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3.44402161 
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-.193652690 
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1.56582916 
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.331719550 









3.1 Muffin-Tin Matrix Elements 



The potential in a muffin-tin at r has an expansion in linear combinations 
of spherical harmonics invariant under that part of the point group leaving 
t invariant: 



V(r) =X^v ht {r T )D ht (V r r r ) 

h 

D ht (r) = ^2a ht (m)Ci hm (r) . 



(18a) 
(18b) 



The utility of referring bases and potentials in muffin-tin spheres to site- 
local coordinates is apparent in (18a). If the site local coordinates of sites 
are constructed so that T> T i = £> T Q _1 for some Q such that Qt = r', then 
the harmonic functions D^t depend only on the symmetry type, rather than 
on each site. The normalization for the spherical harmonic in ( |l8a|) (C = 
\JAtt/ (2^+1)^) is chosen so that Vht{r) is the potential when l) x = 0. 



Combining (|15|) and (18a), the potential matrix is 



{i>i\v\i>i) 



r L 

(XI X n te ( e * > K i)( U u( e i)\ v ht | U u {ej))n u (e 3 , Kj ) 

(L\D ht \L')S L , L] (K 3 ,T~T^k)) . 



h V 
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The matrix element of the Dht is a sum over Gaunt coefficients: 

(L\D ht \L') = ^a ht (m h )g(£', to'; £, to; £ hl m h ) 

rrih 

G(£', to'; £, to; 4, m h ) = / %' m '3^* m C^ mh 



In electronic structure methods using muffin-tin orbitals, the muffin-tin 
energy parameters {eg are usually taken from "^-projected average energies". 
With multiple energy sets, this is a reasonable choice provided that the basis 
set, which uses separate sets, gives rise to bands well separated in energy. The 
^-projected charge, integrated over a muffin-tin sphere, is a sum over cross 
terms between energy sets 



Qi = y^Ql(e»,ej) 



and must be made diagonal in some approximation for the resulting energy- 
and ^-projected energies and charges to be representative. 

Another criterion, particularly useful for states using different sets not 
well separated in energy or for states not having significant occupation is 
to maximize the completeness of the basis. To accomplish this, the energy 
parameter for the low energy state eg(l) can be set to a set of projected 
energy averages, and the energy parameters for the same £ in higher energy 
sets may be chosen so that the radial function has one more node and the 
same logarithmic derivative at the muffin-tin radius, hence 

r 2 dr Mei,r)M«i>r) = > * > 1 • ( 20 ) 



o 

Although this usually generates energy parameters out of the range of oc- 
cupied states (since the logarithmic derivative of semi-core states is usually 
large in magnitude and negative), this choice seems to give a total energy 
close to the minimum with respect to this parameter. This is an example of 
the difference mentioned in the introduction in emphasis between an accurate 
"basis-set" method and a method motivated by a physical model. 

The convergence of the harmonic expansion of the potential in a muffin- 



tin sphere (18a) depends, of course, on the basis, atomic constituents, and 
geometry. Using harmonics through 4 maa . = 6 is usually sufficient, and it has 
never been necessary to go beyond 4 max = 8. 



3.2 Interstitial Matrix Elements 



Overlap and Kinetic Energy: The interstitial overlap matrix can be eas- 
ily obtained from an integral over the interstitial surface (the only non-zero 
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contributions, in a crystal periodic in three dimensions, come from the sur- 
faces of the muffin-tin spheres) and the kinetic energy is proportional to the 
overlap: 

= ( Kj 2 _ K 2)-l^ s 2 /" rffl^^f^.) (21) 
r ^ 

where W(f,g) = fg' — f'g. Basis functions on muffin-tin spheres are given 
in (|l4]), hence 

(ipi\il}j) =J2 s tYl S i,i ( ( K i'" r - T i. fe ) 

r L 

W(Kj(Ki,St),Kt(Ki,S t )) , s , , 



In the limit «j — > /c^, the evaluation of ( |2l| ) requires the derivative with 
respect to k 2 of the structure constant. 



Potential Matrix Elements: The greatest difference between LMTO- 
based full-potential methods is in the way the matrix elements of the potential 
are calculated over the interstitial region. The method being described here 
uses a Fourier representation of basis functions and the interstitial potential 
to calculate these matrix elements. Other approaches for computing these 
elements are described in the literature. |4||| 

A Fourier transform of the basis functions described in Section ^| would 
be too poorly convergent for practical use. However, the evaluation of the in- 
terstitial potential matrix requires only a correct treatment of basis functions 
and potential in the interstitial region. This degree of freedom can be used 
to design "pseudo basis-set" , equal to the true basis in the interstitial region 
although not in the muffin-tin spheres, and have a Fourier transform which 
converges rapidly enough for practical use. We define this pseudo basis set 

by 



&(fc,r) => e^KuiKiAr-Ti-RDi'Y^ir-n-R) (23a) 



R 

JCi{k,t) = /Q(k, r), r > s, s < s T (23b) 



Since rapid Fourier convergence is the criterion for constructing the pseudo- 
basis, it is useful to consider the Fourier integral of a Bloch function with 
wave- number k: 

fo) = - Ve(|fc+ ] |2 „ K2) J^re-^^Wmr) (24) 
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where V c is the unit cell volume. Equation (|2J) is obtained by casting V 2 +k 2 
on the plane wave then doing two partial integrations; surface terms vanish 
due to periodicity. From (|24|) it is evident that the Fourier integral of a 
pseudo-basis satisfying the first criterion (equal to the true basis in the in- 
terstitial region) may be obtained from integral over muffin-tin spheres. If 
in addition, the pseudo-basis is different from a Hankel function only in it's 
parent sphere, the Fourier integral is a finite integral over a single muffin-tin 
sphere. The problem then is to find a function ip such that (V 2 +K 2 )ip has 
a rapidly convergent Fourier integral, vanishes outside a radius less than or 
equal to the parent muffin-tin radius for the basis, and has a value and slope 
equal to K, at this radius. 

A good choice for such a function is obtained by solving 

(V 2 + K 2 )JC e ( K ,r)y L (r) = -<*(-) [l- (-) J y L (r)0(s-r) (25) 

for a radius s < Si ; , and with with eg chosen to match on to K. at s. This is 
easily done analytically. The resulting Fourier transform is 

where N = li+rii+1. The subscript i has been purposely left off N and s 
(see below). 

These coefficients converge like l/g n+4 , provided J^N(\k+g\, s) achieves 
it's large argument behavior, and n can be chosen to optimize convergence. 
Weinert used an analogous construction as tool to solve Poisson's equa- 
tion. He proposed a criterion for the convergence of the Fourier serie ( ^6| ) 
which amounts to choosing the exponent n in Equation ( |26| ) so that \k+g max \s 
would be greater than the position of the first node of Jg+ n +i- We find this 
criterion to be useful provided anisotropy in reciprocal space is accounted for. 
This is accomplished by using the minimum reciprocal lattice vector on the 
surface of maximal reciprocal lattice vectors, rather than simply using g max - 

Notice that this criterion is a criterion for TV = t+n+1. The basis Fourier 
components are simplified, and the amount of information stored reduced, 
by simply using a single argument for all bases; i.e. all bases use the same 
value of N. It is also possible to use a single radius s, less than or equal to 
the smallest muffin-tin radius, since the only requirement is on the pseudo 
bases in the interstitial region. In practice, a few radii are desirable if large 
and small atoms are present in the same calculation, since small radii give 
less convergent Fourier coefficients. In any event, no more than a few radii 
are necessary to handle systems with many atoms. Notice also that local 
coordinates have been left out of (p6f). The resulting potential matrix may 
be easily rotated to local coordinates at the end of the calculation. 

As expressed in (p6|), the Fourier components are products of phases 
e -i(*+s)-T ) which scale like the number of atoms squared (the size of the 
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reciprocal lattice grid grows linearly with the number of atoms), and a func- 
tion of lattice vectors and a few parameters, which scales linearly with the 
number of atoms. The phase factors are simple to calculate by accumulation 
and need not be stored. 

The potential in the interstitial region is similarly obtained from a "pseudo- 
potential" V that equals the true potential in the interstitial region and has 
rapidly converging Fourier coefficients: 



= V(r)\ (27a) 
x I J 



V(r 

V(r) =Y,V(S) D s(r) (27b) 
s 

ges 



The sum in Equation (27b) is over stars S of the reciprocal lattice. 

Integrals over the interstitial region are performed by convoluting the 
potential with an interstitial region step function and integrating over the 
unit cell: 

(^\V\^) X = ($i\V\fa) x = (rP i \e x V\rP j ) c . 

The potential matrix element is calculated by convoluting the convoluted 
potential with a basis, and performing a direct product between convoluted 
and unconvoluted bases. If basis functions are calculated n 3 reciprocal lattice 
vectors, the interstitial potential will be calculated on (2n) 3 vectors. The 
convolution is exact if it is carried out on a lattice containing (An) 3 vectors. 
The size of the set of reciprocal lattice vectors necessary to converge the total 
energy using this treatment of the interstitial region varies from between ~ 
150 - 300 basis plane waves per atom, depending on the smoothness of the 
potential and the convergence required. 

Another way of integrating over the interstitial region, more usual in site- 
centered methods, is to integrate Fourier series over the unit cell and sub- 
tract the muffin-tin contributions with pseudo-bases and pseudo-potential 
expressed as an expansion in spherical waves. The convolution has an advan- 
tage in acting with a single representation, and, given a finite representation 
for bases and potential, the convolution may be done exactly. 

Empty spheres are never used with this scheme. Bases, and the charge den- 
sity and potential are calculated as accurately as necessary using the scheme 
described above and a basis set expanded with tail parameters and energy 
sets has proven to be flexible enough to accurately describe the contribution 
of the electronic states in the interstitial region. 



4 Charge Density 



When a solution to the wave equation at every physical energy is available, the 
charge density may be obtained from a set of energy-dependent coefficients. 
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The spherically symmetric charge density in a muffin-tin sphere, coupled with 
an £ — projected density of states, is an example. In a variational calculation, 
as is being described here, all that is available is a (variational) solution to 
the wave equation at a set of discreet energies, and the charge density must 
be obtained simply from the square of the eigenvectors, or equivalently from 
expectation values of occupation numbers. 

Having calculated a set of eigenvalues and eigenvectors A of the general- 
ized eigenvalue problem, the charge density in the interstitial region is 

n(r) = J2h{S)D s (r) (28a) 
1 s 

h ^ = TrJ2J2 w ^TT f d 3 re-^ r \y2Mk,r)Mnk)\ 2 (28b) 

where Ns is the number of vectors in the reciprocal lattice star S. The square 
of the wave function is obtained by convoluting the Fourier components of tp 
with A, Fourier transforming, and taking the modulus. 
In the muffin-tin spheres the charge density is 

n(r) =Y^n ht (r T )D ht (V T r T ) (29a) 

r ^ <St h 

= EE C/ «'( e ''- r ) MM ( ef ' er ) C/ «( e " r ) ( 29b ) 



el e'l 1 



M ht (e£,e'£') = ^^ ]T a* ht (m h )g(£, m; £', rri; 4, m h ) (29c) 

m h mm' 

^^w nk V T im{e)Vl llm ,{e') 

nk 

Vrim(e) = ^5(e,ei)ntf(e,Ki)St m ,e imi {Ki,T-Ti,k)Ai(nk) (29d) 

i 

The process of calculation is evident in the sequence of equations. 



5 Core States 

Core states, even spherically symmetric complete shells, contribute non-muffin- 
tin components to the interstitial region and to muffin-tin spheres surround- 
ing other sites. Whether it is essential to include this contribution depends on 
the size of the contribution, and any sizable contribution implies that there 
are states being treated as localized which aren't localized within the scope 
of the calculation. Nevertheless, confining states to the core is often useful, 
and including the core contribution to the full potential is not difficult. One 
possibility, the one used in this method, is to fit the part of the core electron 
density to a linear combination of Hankcl functions, and expand this density 
in the interstitial region as a Fourier serie and in the muffin-tin spheres in a 
harmonic series, in the same way the basis functions are treated. 
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6 Potential 

6.1 Coulomb Potential 

The Coulomb potential is obtained by first calculating the Coulomb potential 
in the interstitial region, then, using the value of the interstitial potential on 
the muffin-tin sphere, calculating the potential in the spheres by a numerical 
Coulomb integral of the muffin-tin electron density for each harmonic. 

The interstitial Coulomb potential is calculated in a way similar to that 
suggested by Weinert |12j . Express the electron density as 

n(r) = h(r) + /J(n(r) — n(r))0(s t — r T ) (30) 

Rt 

where n is the squared modulus of the pseudo-eigenvectors, which is equal to 
the true electron density in the interstitial region. The hrst term on the right- 
hand side of ( |30| ) has, by construction, a convergent Fourier series. The second 
term is confined to muffin-tin spheres. To calculate the Coulomb potential in 
the interstitial region, this term may be replaced by any density also confined 
to the muffin-tin spheres and having the same multipole moments. If a charge 
density satisfies these requirements and also has a convergent Fourier series, 
the Coulomb potential in the interstitial region may be easily calculated from 
the combined Fourier series. Such a charge density can be constructed in a 
similar way to that detailed for the pseudo-bases. Construct a pseudo charge- 
density satisfying 



n(p)(r) ^^AW^,^)^^^) (31a) 

Rt h 

< P t(r) = c w (f ) ih (l (f ) 2 T0(s t -r) (31b) 

\ St ' > > St ' ' 

= J d 3 rr e T D* ht (V T r T )(fi ip) {r)-n(r)+r~i(r)) . (31c) 

This charge density has Fourier components 

^ (r) = E E e-^(-iY»D ht (V Tg) ^ (Q^;>^« W) 
h 

(2(4+n+l) + l)!! 



h 



9 th Ji h+n +x{9,s t ) (32) 
where the multipole moments Q are defined by 

Qht{n} = ?^±I / 4»D ht (r T )n{r) d\ T (33) 

4?T J St> r T 

The Fourier components fv- p '(r) converge like l/g n+2 provided j^„-\-i at- 
tains it's asymptotic form. The exponent n is chosen using the same consid- 
erations as for the pseudo-basis set. 
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The Coulomb potential in the interstitial region is then given by 

V c (r) =V c (r)\ 
i \x 



E 



r 



(34) 



From the Coulomb potential in the interstitial region follows the Coulomb 
potential on the surface of the muffin-tin spheres. The coulomb Potential 
inside the muffin-tin spheres is 



V {c \r) 



r T <s t 



Y,D ht {V T r T ) 



s t r th 4 7rr ' 2 7l/j ( r ) f 



th+i 2£ h +l 



-dr' 



(35) 



ATrr' ih+2 n h {r') a fr\^ 



24+1 



(j) 



where 



24+1 

47T 



drD* ht (V T r)V^(r) 



(36) 



is the harmonic component of the potential on a sphere boundary. 



6.2 Density Gradients 

Gradients of the electron density are needed for the evaluation of gradient 
corrected density functionals. These functionals depend on invariants (with 
respect to the point group) constructed from density gradients {e.g. |Vn| 2 ). 
This reduces computation significantly in the muffin-tin spheres, for if / and 
g are invariant functions (i.e. f(r) = J2h fh{r)Dh(r)), and d = V/ • Vg, 
then d(r) = J2h d h {r)D h (r) with 

P^Mr) = J2 E fi k \r)g%'\r)I(kk';hh') (37) 
h h,h' k,k'=±l 

where the set of parameters / is easily calculable from 3j and 6j coefficients 
and integrals over the harmonic functions Dh, and 

/(*)_ 47r frf'-th k = l 

h ~ 21+1 \rf +4 + 1 k = -l [ ' 

and similarly for g. 

Gradients of the interstitial charge density, represented as a Fourier se- 
ries, are poorly represented by differentiating the series term by term. A 
stable representation of the density gradient that converges well is obtained 
by defining the derivative as the difference between adjacent grid points, di- 
vided by twice the grid spacing as suggested by Lanczos. [jl3| This is equivalent 
to differentiating, term by term, the Lanczos-damped series for the charge 
density. 
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7 All-Electron Force Calculations 
7.1 Symmetry 

The set of internal forces acting on the atomic sites of a crystal is a symmetric, 
discrete function of atomic coordinates and has a spherical expansion on the 
crystal sites with the same coefficients as continuous symmetric functions 
( js"a|) and (|8b|). Since forces are vectors, their representation has 1=1, and if 
a site has no invariant harmonics with I = 1, there is no force on that site. 
So the force on an atomic site may be expressed as 

f(r)= £ fhtJ2 

U T (39) 

h:l h = l m 



where the coefficients a are as in (8b), the e m are spherical unit vectors, 



[[L4| and U T is the transformation to local coordinates for spherical vectors. A 
force calculation is, as much as possible, a calculation of the set {fht}', The 
size of this set is often much smaller than three times the number of atoms. 
The displacements of atoms allowed by symmetry also have the form of (p9|) : 



St = £ 5r ht £ a m e m U T (40) 



h-i h =i 



Minimizing the energy with respect to the atomic positions is a process 
of finding the set {Stm} that gives fht = 0. 



7.2 Force Calculations 

The calculation of forces in an all-electron method has been nicely described 
by Yu et al. for the LAPW method. In addition to the terms discussed in 
that paper, a force calculation using a site-centered basis has the additional, 
and significant, complication that the bases depend on atomic position not 
only through augmentation but also through parentage. 

The contributions to the total force on a site in an all-electron calculation 
follow directly from a derivative of the LDA total energy with respect to 
atomic positions. The terms listed by Yu et al. are 1) a "Helmann-Feynman" 
term, dE/dr, which accounts for the explicit dependence of the energy func- 
tional on atomic positions, 2) an "Incomplete Basis Set" (IBS) term, which 
arises when derivatives of basis functions aren't contained in the space cov- 
ered by the basis set, 3) a core-correction term, arising because core states 
are calculated using only the spherical average of the potential, and 4) a 
muffin-tin term, a surface term arising from the change in integration bound- 
aries when atoms are moved and the discontinuity of the second derivative 
of basis functions across muffin-tin boundaries. There are two other terms 
to consider. The first arises when a calculation isn't fully self-consistent, and 
has the form — J y (V ou t — Vi n )dn(r)/d,T, where V ou t and VJ n are output and 




-0.04 



0.9 



1.0 

V/ V 



1.1 



exp 



Fig. 1. The deviation of the internal coordinates of rhombohedral BaTi03 from 
ideal, calculated using all-electron force calculations as a function of volume with 
both LDA (open symbols) and GGA (filled symbols) exchange-correlation func- 
tions. The grey filled symbols are experimental points jl6). The LDA equilibrium 
volume is .958 V exp ; the GGA volume is 1.037 ~V exp - The energy was also minimized 
with respect to the rhombohedral angle at each volume. 



input potentials. The second term arises from the way in which Brillouin zone 
integrals are done. Whether by quadrature or linear interpolation, the result 
is a set of weights (occupations) multiplying quantities evaluated at discrete 
Brillouin zone points. The terms listed above do not take into account the 
change of weights with atomic positions. 

The evaluation of the IBS term in a method using site-centered bases is 
significantly more involved than in the LAPW method. This term has the 
form 

f ibs = -^2w nk ^2A* nk ( (ipi\H - e nk \di()j/dT) 

nk ij 

+ (<tyi/d,T\H - e nk \ipj))Aj tn k (41) 

where the A are eigenvectors. Both LAPW and LMTO methods have a de- 
pendence on atomic positions through augmentation (the expansion of the 
basis set in atomic-like spherical waves) in the muffin-tin spheres, and both 
methods have an implicit dependence of basis functions on atomic positions 
through self-consistency, a term largely ignored and usually negligible. A site- 
centered basis, however, depends on atomic: positions also through it's parent 
site (the site it's centered on). The contribution from augmentation is fairly 
easily accounted for at the density stage of a calculation, after integrals over 
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the Brillouin zone have been done. The parent contribution, however, requires 
evaluation at the part of the calculation where eigenvalues and vectors are 
obtained, which makes its calculation time consuming. 
There are four types of contributions to d^jdr: 

-^ i (k,r)=i(j5W+S? ) +8? ) +6W)i> i (k,r) (42) 
tf( 1 Vi(*,r) = 0(r e l)5{n,T) ptpi(k,r) (43) 
8^ipi(k, r) =5(Ti,T)^2 @(st>-r T ')U t >L(ei, r r ,)n vt {e u m) 



X (-iV T S i , ii (/t i) T / -T 4s fc)) (M) 

S^ipi{k,r) = 0(s t -r T )^2p\JtL{e l ,r T ){2u(ei,K l )S L . Ll (K i ,T-T i ,k) (45) 



6^ipi{k,r) = -Q(s t -r T ) ^ Utz,(e», r T )Q ti (ej, m) 

L 





-iV T B LtLi (Ki,T-Ti,k) 



(46) 



where p is the momentum operator — iV. The first two terms, Equations 
( (43|) and (HJ), are parent terms, changes in a basis due to a change in the 
site the basis is centered on. The first term, Equation ( ^3|), is the derivative 
of the wave function in the interstitial region (Equation (|l0|) with respect to 
its parent site. Since the gradient of a solution to the Hclmholtz equation is 
a solution to the Helmholtz equation, matrix elements (^i\pij)j} x and — 
V 2 \ptpj ) I are calculated as integrals over the surface of the muffin-tin spheres. 
As in Equation (^||), when interstitial region tail parameters are the same, 
the evaluation requires k 2 derivatives of structure functions. Working out this 
contribution proceeds as in Equation (S), although arriving at a finite form 
requires identities such as 

Tb,k)G(£b-l, rn b -fi; 4, m b ; 1, /i) 

ft 

- B i a m a ,i i +im b -n( K b,T a -T b ,k)g(£ b +l ) m b -{i; i b , m b ; 1, fi 
T a ( B la+lma+lllbmb (K b ,T a ~T b ,k)g(l a , m a ; 4+1, m +/x; 1, fi 

B e a -i ma +n,e b m b ( K b,T a -T b ,k)g(£ a , m a ; 4-1, m a +fi; 1, ^ 



ft 




As 




Fig. 2. Relaxation of a silicon 65 atom supercell containing a vacancy, a Si inter- 
stitial, and an As interstitial. Of the 106 internal coordinates in this cell, 104 were 
allowed to relax (2 coordinates were fixed to fix the center of mass of the crystal). 
The calculation used a simple Broyden's method to zero atomic forces. 



Potential matrix elements are calculated using Fourier series as in 

Sect. 3.2 with gradients taken as discussed after equation (|38|). 

The second term, equation (|4^), is the analog of the first term in the 
muffin-tin spheres; i.e., this term is the derivative of a basis with respect 
to its parent site evaluated in the muffin-tin spheres. This term requires the 
gradient with respect to atomic positions of the structure function B. This 
gradient is easily obtained from the structure function itself: 

8 



U — T — T 1 



= ^ie^ r sg /m ,( K ,T-r',fe) 

A" 

-K 2 g(£, m; £-1, m+n; 1, fi)B e _ lm+ ^ A , m ,(K 7 T- 
t-t'/O 
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If convergence with respect to I on the left hand side of the structure func- 
tion is sufficient for the energy, terms in i max + 1 in Equation (^8|) may be 
neglected in evaluating forces. As stated above, the evaluation of these terms 
is somewhat time consuming. 

Examples of the use of forces for structural relaxation are given in Figures 
and |^. Figure [l] shows deviations from ideal lattice positions calculated for 
rhombohedral BaTi03 as a function of volume compared to experiment ]l6| . 
The rhombohedral angle was also relaxed at each volume in this calculation. 
The Ti coordinate is a displacement along [111]. The oxygen displacements 
Ax are along face diagonals while Az is toward the cell center. These calcu- 
lations included Ti 3s and 3p and Ba 5s and 5p along with the usual valence 
bases in a single, fully hybridizing basis. At convergence, forces on internal 
coordinates were less than 1 mRy/Bohr. Figure ^ is a calculation of struc- 
tural relaxation of As- vacancy-interstitial complex in Si. To a sixty-four atom 
Si supercell was added an As impurity at a tetrahedral interstitial position 
and a Si interstitial at an exchange position both surrounding a vacancy. The 
crystal, far from equilibrium, was then allowed to relax. Two internal coor- 
dinates (of a total of 106) were fixed to fix the center of mass of the crystal. 
The energy was minimized with respect to the other 104 internal coordinates 
by zeroing the forces (to with 1 mRy/Bohr). The forces were zeroed using a 
simple Broyden's method. 



8 Conclusion 

In this article we have described our highly accurate full-potential LMTO 
method for solving the Kohn-Sham equations. In particular, we have shown 
that by dividing the crystal space into non-overlapping "muffin-tin" spheres 
and an interstitial region, we can compute the charge density or the poten- 
tial without any shape approximation, thus eliminating any need for empty 
spheres which are necessary in other LMTO implementations when the crys- 
tal is not closely packed. Another feature of our implementation is that we 
can describe multiple principle quantum numbers within a single, fully hy- 
bridized basis set. This is accomplished simply by using functions <f> and <fi 
calculated with energies {e n i} corresponding to different principal quantum 
numbers n to describe the radial dependence of a basis in the muffin-tin re- 
gion. In the interstitial region our method uses "multiple k" basis sets, for 
a better description of the interstitial charge density. Highly accurate charge 
density can be obtained by systematically increasing the number of varia- 
tional parameters k for each angular momentum of the basis set. 

The potential in a muffin-tin sphere at r has an expansion in linear com- 
binations of spherical harmonics invariant under that part of the point group 
that leaves atomic positions invariant. The evaluation of the interstitial po- 
tential matrix only requires a correct treatment of basis functions (and poten- 
tial) in the interstitial region. We have used this degree of freedom to design 
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"pseudo-basis functions" , equal to the true basis functions in the interstitial 
region and are smooth functions in the muffin-tin region, with the require- 
ment that their Fourier transforms converge rapidly enough for practical use. 

The set of internal forces acting on the atomic sites of a crystal is a sym- 
metric, discrete function of atom coordinates and has a spherical expansion 
on the crystal sites with the same coefficients as continuous symmetric func- 
tions. The total force on a site is given by the derivative of the LDA total 
energy with respect to the atomic position. Our implementation of the forces 
is in many ways similar to that of Yu et al. for the LAPW method fl5|| . 
Because our basis set is a site-centered one, we are required to compute addi- 
tional terms, which can be time consuming. These contributions to the forces 
are non existant in plane-wave based methods, such as the pseudo-potential 
method. In addition to the "Helmann-Feynman" term, which accounts for the 
explicit dependence of the energy functional A on atomic positions, the other 
contributions are: (1) an "Incomplete Basis Set" term, (2) a core-correction 
term, (3) a surface term arising from the change in integration boundaries 
when atoms are moved, (4) a term which arises when the calculation isn't fully 
self-consistent, and (5) a term arising from the way in which the Brillouin 
zone integrals are performed. We have showed that the forces are accurate 
enough to relax atomic structures. As examples, forces have been used to 
optimize the internal coordinates of rhombohedral BaTi03 as a function of 
volume and the geometry of a 65 atom As, vacancy, and interstitial defected 
Si supercell. Where experimental results are available, good agreement is 
obtained. 
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